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ABSTRACT 


Littoral waters when compared to the open ocean create an environment of greater 
reverberation with acoustic energy scattering from the sea surface, bottom, topographic 
features, and regions that lack homogeneity within the volume. If the ocean surface is rough 
on the scale of an acoustic wavelength, considerable scattering can occur that can 
significantly influence coherent propagation. Because the rough surface is also evolving 
dynamically, such scattering can introduce Doppler shifting and spreading of the acoustic 
pulse spectrum. This thesis builds upon prior efforts in ocean acoustic modeling and is 
focused on examining surface scattering and its affect upon coherent propagation. The 
dynamics/physics associated with sea surface scattering are explored in detail and 
mathematical relationships are developed and employed in revisions to the Monterey Miami 
Parabolic Equation (MMPE) model. The thesis provides background information associated 
with the MMPE and highlights earlier work related to surface scattering. It presents a formal 
analysis of an exact surface scattering approach in the context of a continuous wave (CW) 
benchmark exercise and the Doppler shifts associated with a dynamic rough surface. It 
expands on prior rough sea surface work to include modeling based on an empirical fetch- 
limited ocean wave spectrum and compares modeling results with measured data. Interest in 
broadband pulse propagation in shallow water is increasing with the need for improved 
active sonar systems and with the growth of applications such as underwater acoustic 
communications. 
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I. INTRODUCTION 


A. THESIS OBJECTIVES AND MOTIVATION 

The principal objectives of this thesis are to demonstrate the capability to model 
the complex physical interactions associated with acoustic broadband pulse propagation 
reflecting from a rough sea surface. The sea surfaces modeled include cases for static 
and dynamic simple surfaces and more complicated realizations of static surfaces. 

Littoral waters when compared to the open ocean create an environment of greater 
reverberation with acoustic energy scattering from the sea surface, bottom, topographic 
features, and regions lacking homogeneity within the volume. If the ocean surface is 
rough on the scale of an acoustic wavelength, considerable scattering can occur that can 
significantly influence coherent propagation. Because the rough surface is also evolving 
dynamically, such scattering can introduce Doppler shifting and spreading of the acoustic 
pulse spectrum. Interest in broadband pulse propagation in shallow water is increasing 
with the need for improved active sonar systems and with the growth of applications such 
as underwater acoustic communications. Also, determining sea surface roughness 
presents an inverse problem of interest in which acoustic waves can possibly be used as a 
tool to understand the physical environment. 

B. SUMMARY 

This thesis builds upon prior efforts in ocean acoustics modeling and is focused 
on examining surface roughness and its affect upon coherent propagation. The 
dynamics/physics associated with ocean sea surface roughness are explored in detail and 
mathematical relationships are developed and employed in revisions to the Monterey 
Miami Parabolic Equation (MMPE) model. 
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Section II covers background information in detail associated with the MMPE and 
highlights earlier work related to surface scattering. Much of this earlier work centers on 
the use and upgrade of acoustic propagation codes based on the split-step Fourier 
parabolic equation algorithm introduced by Hardin and Tappert (1973). As an upgrade, 
Tapper! and Nghiem-Phu (1985) developed an algorithm to compute exact surface 
forward scatter when the surface interface is not assumed to be flat. An approximate 
rough surface scatter was later developed and published in a University of Miami 
Parabolic Equation (UMPE) technical report (Smith and Tappert, 1993). Other published 
references to the approximate scatter approach are found in two theses (Wei Ei, 1993, and 
Jun He, 1997). 

Section II also describes implementation of the MMPE model and details such as 
operation of the model can be found on the Ocean Acoustics Eibrary website at 
http;//oalib.saic.com/PE/index.html. Section II concludes with a literature search of 
earlier work showing how surface waves act as moving diffraction gratings, scattering 
low frequency sound in selected, well defined directions and imparting a unique Doppler 
shift in each direction. 

Sections III and IV cover a formal analysis of an exact surface scattering 
approach in the context of a continuous wave (CW) benchmark exercise and the Doppler 
shifts associated with a dynamic rough surface. Section III implements into the MMPE 
model a static rough surface forward scatter subroutine similar in function to what has 
been developed for the UMPE model. The influence of a simple sinusoidal surface 
perturbation is examined and Doppler effects from a dynamic ocean surface are 
computed. Section IV sets up the model to calculate Doppler shifts in the presence of a 
dynamically moving surface and verifies that the updated model performs adequately. A 
benchmark exercise is conducted with the results showing the expected frequency spectra 
at different depths and different ranges. 

Section V expands on the rough sea surface work to include modeling based on an 
empirical fetch-limited ocean wave spectrum and compares modeling results with 
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measured data. The Pierson-Moskowitz spectrum is presented showing energy spectra 
for fully developed waves at various wind speeds. Next, an experiment conducted in 
1997 (HFA97) measuring cause and effect between the ocean environment and acoustic 
propagation is analyzed. The MMPE model is updated with the JONSWAP wave 
spectrum with results clearly showing surface reflected paths dispersing acoustic energy 
at increasing rates as wind speed increases, and also clearly showing how the direct path 
(non-surface interaction) remains unchanged. MMPE Model statistics are computed for 
standard deviation of arrival angle and arrival time and compared with both HEA97 
measured data, and data from another model (BEEEHOP/JONSWAP). 

Conclusions and recommendations are discussed in Section VI. The results 
presented show that surface scatter can significantly affect arrival time, arrival angle, and 
frequency spread that leads to degradation of acoustic signal coherence. Some 
recommendations for future work include examining Doppler shifts for the dynamically 
evolving rough surface case, improving upon data processing MATEAB Algorithms, and 
validating algorithms with additional measured data. 
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II. PREVIOUS WORK 


A. BACKGROUND AND DESCRIPTION OF MMPE MODEL 

I. Analytical Development of MMPE 

Hardin and Tappert (1973) are two of the earlier contributors to have applied a 
parabolic approximation of the acoustic wave equation to predict underwater sound 
propagation. Since then, parabolic equation models have become the models of choice in 
cases where the environment varies with range. 

The predecessor to the Monterey Miami Parabolic Equation (MMPE) model 
(Smith, 2001) is the University of Miami Parabolic Equation (UMPE) model developed 
as a research model under the guidance of Professor Ered Tappert (Smith and Tappert, 
1993). The MMPE is based on most of the same approximations and numerical 
algorithms as the previous UMPE model with a few exceptions. The most notable change 
includes the centered-step scheme of the split-step Eourier algorithm. This change 
improves accuracy in the range-step calculation by an order of magnitude while adding a 
fractional amount to the overall run time. 

The MMPE model is derived by beginning with the definition of the Helmholtz 
wave equation in cylindrical coordinates, 

V^p{r,z) + Kln^{r,z) p(r,z) = 0 . (2.1) 

Equation (2.1) is then factored by introducing the operator notation 

1 

Q„^={p + s + iy , ( 2 . 2 ) 


where 


5 



(2.3) 




1 a 

dz^ ’ 


s = n^ , 


and 



(2.4) 


(2.5) 


The outgoing pressure field may then be defined in terms of the parabolic 
equation (PE) field function, y/, according to the following equation, 


p{r,z) 



Qop z) e" 


( 2 . 6 ) 


where y/ satisfies a parabolic equation of the type 


dy/ 

dr 


-ik^y/ + ik^Q^^y/ = -ik^H^^y/ . 


(2.7) 


The split-step Fourier algorithm (Hardin and Tappert, 1973) is then applied by 
separating the “Hamiltonian” operator, H , into the sum of “kinetic energy” and 

“potential energy” operators, and U , such that the field is marched outward in 
range according to 


^(r + Ar,z) =e 




(r+Ar,z) 


f 


FFT\ 


-ik„ArT(kJ 


IFFT 


-ik„ — U„ 
° 2 


(r,z) 


¥ir,z) 


( 2 . 8 ) 


In the MMPE model, the expressions for the operator approximations employ the wide 
angle PE (WAPE) forms (Thompson and Chapman, 1983), 
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(2.9) 


U „p{r,z)^-[n{r,z)-\] 


and Top{k ) = 1- 



(k ^ 

2 

1- 







Finally, in terms of the surface reflection, the pressure release boundary condition 
is imposed at each range step by forcing the odd symmetry constraint 

y/{-z) = -y/{z) , (2.10) 

which requires the use of an “image ocean” solution when using the full FFT form (rather 
than a simpler sine transform for z>0). Further details of the MMPE implementation can 
be found in the review article by Smith (2001). 

In the case of a rough surface, the pressure release boundary condition must still 
be satisfied, but now it is imposed at the position of the displaced surface defined by 

z-; 7 (r)= 0 , (2.11) 

such that 

y/iz = ri{r) = 0. (2.12) 

Tapper! and Nghiem-Phu (1985) showed that this can be achieved by defining the 
Uop operator (dependent on the environmental index of refraction) as an even function 

and the field function y/ as an odd function about the displaced surface interface, and 
then solving two parabolic equations, one for the real ocean and one for the image, 
defined by 


real ocean: 


dy/ 

dr 


i^a^op+Uop)y/, 


Z> ri{r). 


(2.13) 
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and 


image ocean: 


dy/ ^^drj dy/ 
dr dr dz 




z > rjir). 


(2.14) 


It can be shown that 


y/{r-z + 2ri{r)}=-y/{r, z), 

(2.15) 

y/(r,7j(r))=0, 

(2.16) 


as required. 


The two parabolic equations are transformed into the single form previously 
defined by introducing a new field function extending over both real and image ocean 
depths defined by 


y/{r,z)= 


¥{r,z). 


/ rs / W 

-e i//[r-z-\-2?j[r )), 


Z>ri{r) 

Z<Jl{r) 


(2.17) 


This new field function then satisfies the previous parabolic equation form at all depths, 
but with an altered U , i.e. 


dy/ 

dr 


f 


ikr. 




Zp+Uop 

V y 


w . 


(2.18) 


where 


U op{r,z)= 


\Uop{r,z), 

Ujr,z)-2{z-rj) 


d^rj 


z>ri{r) 

z<ri{r) 


(2.19) 
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2. Implementation of MMPE 


The code for the MMPE was developed in Eortran and is accessible with a Eortran 
compiler such as Digital Visual Eortran Development Studio. Details as to the makeup 
and operation of the model can be found on the Ocean Acoustics Eibrary website at 
http ;//oalib. saic. com/PE/indel. html. 

The model’s executable file, mmpe2dbb.exe, is compiled from six source files 
and calls seven ascii input files during operation. These files contain information on run 
parameters, source parameters, and all necessary environmental information. 

When the program is run, a single output binary file is created in which is 
contained a header and the vast majority of the remaining file is the PE field function 
y/i^r, z ). The mmpe2dbb executable is a broadband version of the 2-D PE model. 

Two post-processing files developed in MATEAB code process the binary file 
data. To initialize the output processing one program (peoutl) is run (in MATEAB) 
which reads the header and provides file identification for further analysis. Next, a 
second program (peout2) is run (in MATEAB) which provides menus for selecting 
various options for processing the data. Eigure 2.1 shows the MATEAB display results of 
computing data for a single radial with a source frequency of 400 Hz, a source depth of 
50 meters, a water depth of 250 meters, and at a range of 10 km for a simple, range- 
independent environment. 

B. LITERATURE SEARCH 

Warfield (1981) describes Doppler shift scattered from the ocean surface well in 
his article titled “Doppler Shifting of Surface-Scattered Reverberation”. The underlying 
physical mechanism is that gravity waves on the ocean surface act as moving diffraction 
gratings, scattering low-frequency sound in selected, well-defined directions and 
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Figure 2.1 Transmission Loss vs. Depth and Range 


imparting a unique Doppler shift in each direction. One aspect of the reverberation 
problem that appears to be well in hand is the Doppler shift of CW plane waves scattered 
from the ocean surface and sensed by a directional receiver. 

Liebermann (1963) showed that backscatter of sound waves in air from a moving 
water surface preferentially selects a wavelength from the water surface spectrum 
according to the familiar diffraction-grating equation. Since each wavelength has a 
characteristic speed, the backscatter has a characteristic Doppler shift (even though no 
basic wave physics analysis was presented to indicate why the diffraction-grating 
equation should work); his experiment was probably the breakthrough result on the 
problem. 

Marsh (1963) presented an analysis shortly after Liebermann of the Doppler shift 
of boundary reverberation based on scattering theory. His analysis gave an expression for 
bistatic reverberation spectra in which the Doppler shift is determined according to a 
diffraction-grating relationship. 
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III. UPDATED MMPE EOR SUREACE SCATTERING 


A. STATIC ROUGH SURFACE MODELING 

The first modification to the MMPE code was to provide a representation of a 
static rough surface. A rough surface forward scatter subroutine (ZSGEN) similar in 
function to what has been developed for the UMPE model was added to the MMPE to 
calculate a roughness spectrum. Eor the purpose of computing the acoustic field in two 
dimensions (depth and range), only a 1-D roughness spectrum along the track of interest 
was needed. 

The subroutine ENVPROP was modified to perform a calculation of an image 
ocean potential function. The output data from ZSGEN provides the input for this 
calculation. Again, the modifications are similar to what has been implemented 
previously in the UMPE. The surface is treated as a perfect reflector due to a pressure 
release boundary. With this method, we assume an identical image ocean overlays the 
real ocean for negative values of depth and, furthermore, the acoustic field is exactly 
equal but of opposite sign in the image ocean. 

The main routine, PEMP, was updated to compute exact rough surface scatter and 
transform the results to the physical space domain. When the model is run, the user is 
able to explicitly define the rms surface roughness (see Pig. 3.1). 

The model was run with an rms roughness of 1, 5, and 10 meters. The input 
parameters are as follows: source depth of 50 m, center frequency of 400 Hz, maximum 
range of 10 km, maximum depth of 400 m, and single radial output. 
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Figure 3.1 Model Running Display 


Figures 3.2, 3.3, and 3.4 show results for an rms roughness of 1, 5, and 10 m 
respectively. Note the increase in near surface attenuation as the rms roughness increases. 
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Figure 3.2 Transmission Loss for rms Roughness of 1 m at 400 Hz 
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Transmission Loss (dB re 1m) with rms roughness of 5 meters 
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Figure 3.3 Transmission Loss for rms Roughness of 5 m at 400 Hz 


Transmission Loss (dB re 1m) with rms roughness of 10 meters 


0 2 4 6 8 10 

Range (km) 


Figure 3.4 Transmission Loss for rms Roughness of 10 m at 400 Hz 
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B. SCATTERING FROM A SIMPLE SINUSOIDAL SURFACE 


As a test of the exact forward scatter model, and to introduce the ability to 
compute Doppler effects from a dynamic ocean surface, we examine the influence of a 
simple sinusoidal surface perturbation. A sinusoidal perturbation will create Bragg 
resonant (diffraction grating) scattering, where the reflected field will be concentrated 
around the Bragg lines (see Fig 3.5) which satisfy the following: 

cos6*^ „ = cos6',. ± n, (3.1) 

where Z = acoustic wavelength, = surface wave wavelength, 0. = incident angle, and 
0^^ s reflected angles. 



n =-2 

Figure 3.5 Depiction of Bragg scatter 

For the test case conducted, we set / = 1500 Hz, c=1500 m/s, Z = l m, and 
= 20 m. 


Then 




(3.2) 
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If 6*, = 30°, the Bragg lines occur at 
^r.O = 30 ° 

6', 1 = 23.6° = 35.3° 

0^2 = 15.0° 6», _2 = 40.0° 

^r ,3 ~ complex - 44.3° 

It was determined empirically that the updated MMPE needs good depth 
resolution of the surface displacement (Z). For a maximum displacement of Z, the depth 
mesh should sample at least by Z/5 (or better). For small surface roughness this puts a 
large computational burden on the model because it requires very fine depth increments. 
Further, it hasn’t been determined what minimum value of roughness can be ignored 
(though it should be dependent on the acoustic wavelength). Also, the model was written 
in such a way that the surface displacement is simply rounded to the nearest depth mesh. 

To observe Bragg scattering requires an incident plane wave from which the 
reflected/scattered energy is observed propagating in specific directions. This was 
accomplished by assuming the source was a long (1000 m) vertical array steered towards 
the surface (D/E of - 30°). 

The subroutine ENVPROP was modified to input from the operator the surface 
wave speed and the time of surface realization. Sinusoidal surface displacement was then 
calculated. The model was run first simulating a flat surface (0 m roughness). It was then 
run with 5m rms roughness. The MATFAB output subroutine PEOUT2 was also 
modified to display the upper part of the water column only (the better to observe Bragg 
scattering). 

The sound speed profile file was set to isothermal conditions at 1500 m/s. The 
model was run with settings that limit the output of depth data to 400 meters in order to 
save on processing and storage time. The maximum range was limited to 4 km. 
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Figure 3.6 shows results from the interaction in the upper part of the water 
column of the plane wave reflecting from a flat (0 m roughness) surface. Figure 3.7 
shows modeling results using a 5 m amplitude sine wave surface roughness with a 
surface wavelength of 20 meters. The absolute value of the difference between the 
complex pressures computed in each case is presented in Fig. 3.8. Figure 3.9 is a plot 
with the same data as Fig. 3.8 but now the Bragg lines have been superimposed on the 
figure at the proper Bragg angles. 


TL (dB re 1 m) with Flat Surface 



Range (km) 

Figure 3.6 Flat Surface Interaction from Upper Part of Water Column 
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TL (dB re 1m) with 5m Sinusoidal Surface 
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Figure 3.7 5m Roughness 
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Figure 3.8 Absolute Value Difference 
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ABS(Press(5m) - Press(Om)) w/Bragg Lines ^ 
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Figure 3.9 Bragg Lines Superimposed on Difference Plot 
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IV. SIMPLE TEST OE DYNAMIC SUREACE-INDUCED 

DOPPLER SPREAD 


The focus in this section is setting up the model to calculate Doppler shifts in the 
presence of a dynamically moving surface and to verify that the updated model performs 
adequately. 

A. INCORPORATING SIMPLE, DYNAMIC SINUSOIDAL SURFACE 


For a moving surface wave with frequency co^, the surface wave number is 


k 


w 



and the surface displacement phase speed will be 


V 


w 


'In 


(4.1) 


(4.2) 


The scattered beams will then undergo a Doppler shift according to (Medwin and Clay, 
1998) 


(O^ = (O, 


1 + — (cos 0^ - cos 6*,) 
c 


= (O; 


c X... 


= (O; 


c In 


(4.3) 
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or with / = —, 
In 


fr =/,■ 


f 

l±n^ 

/, 


= fi ±«/w > 


(4.4) 


so 




2n 2n 



(4.5) 


Recall the following for the test case in Chapter III: /=1500 Hz, c=1500 m/s, 
A = lm, and = 20m. To make this noticeable at f. = 1500Hz, we desire «50 
Hz. This requires that = 1000 m/s! This is not realistic, but will be used to test 

the model. 

To ensure adequate sampling of a 20 m wavelength moving surface wave, we 
should sample roughly every 2 - 4 m. Since the surface wave is moving at 1000 m/s, this 
suggests a minimum sampling rate of = 250 Hz. Furthermore, since the Doppler is 

» 50 Hz, we should sample in frequency at least every 5-10 Hz, thereby requiring a time 
sample length of at least T =0.256sec. 


Based on these requirements, the model was run 128 times with time steps of 
At = 2 msec (/^ =500Hz) for a maximum time span of T = 0.256sec. At each time 
sample, the surface displacement was computed from 


//(r) = Asin 




(4.6) 
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where A is the surface wave amplitude, set to 5 m. The resulting 128 point time series 
data were extracted at various positions within the field and Fourier transformed to 
observe Doppler shifts. 

Next, data was collected with the model running 32 times and incrementing in 
steps of 0.003125 seconds from t = 0 to t = 0.\ seconds. The rms roughness was kept 
constant at 5 meters. The surface wave speed was kept constant at 1000 m/s. Figures 4.1 
through Figs. 4.5 represent trial numbers 1, 8 , 16, 24, and 32 with t=0, t=0.02581, 
t=0.048387, t=0.074194, and t=0.1 seconds respectively. Observe the Bragg scatter by 
noting how the changes in the reflected rays (as time progresses) travel in distinct 
directions. 
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Figure 4.1 Trial #1, t=0 


21 








Depti (m) 
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Figure 4.2 Trial #8, t=0.022581 
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Figure 4.3 Trial #16, t=0.048387 
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Figure 4.4 Trial #24, t=0.074194 
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Figure 4.5 Trial #32, t=0.1 
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B. REALIZATION 


Data was extracted at single ranges of 2.5 and 3.0 km. This produced pressure 
matrices in the MATLAB workspace as functions of depth and frequency. From these 
pressure matrices, data vectors were extracted at depths of 100, 150, 200, 250, 300, and 
350 m. An inverse FFT followed by an FFT shift was performed on the data vectors and 
finally the absolute value was computed prior to displaying the results. 

To automate the process a shell program was developed in MATLAB that cycles 
128 times calling the MMPE update program and saving the wave speed and elapsed 
time, t, each cycle for use by the environmental propagator subroutine, ENVPROPl. The 
time t is incremented with the formula, t=(N-l) * Dt, where Dt is set to 0.002 seconds. 
The calculation for the surface displacement. 


ri{r) = A 


sm 


^{r + vj) 


(4.7) 


is computed by ENVPROPl. When the program is run (approximately 8 hours at 2.2 G 
Hz processing speed) 128 bin files are generated. 

Next, PEOUTl and PEOUT2 were combined into one MATEAB program, 
PEOUTllEN, and treated as a function. A main program, PEOUMAIN, was developed 
which loops 128 times performing the computations as described previously and outputs 
the final plots. Eigs 4.6 through 4.17 show the results with the frequency spectrum 
varying from + and - 150 Hz in 50 Hz steps. 
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c. 


EXAMINING RESULTS 


According to the Bragg scatter conditions, there should be Doppler signatures at 
±50Hz, ±100Hz, and potentially even a contribution at ±150Hz. The results in Figs 4.6 
through 4.17 show energy predominately at ±50 Hz, and ±100//z, with a small amount 
of energy occurring at bins ±150//zin Figs 4.6, 4.8, 4.10, 4.12. 4.14, and 4.16 (all at a 
range of 2.5 km). Thus, the updated model has been verified effective at computing 
Doppler spread due to surface motion by using the exaggerated surface wave speed of 
lOOOm/s. 


X 10^ Figure 1. Ocplh=100 m. Range=2.5 km 



Figure 4.6 Frequency Spectra at a depth of lOOm and range of 2.5 Km 
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Figure 2. 0epth=100 m. Ran9e=3.0 km 


X lO* 



Figure 4.7 Frequency Spectra at a depth of 100m and range of 3.0 Km 


Figure 3. Dep1h=150 m. Ran9e=2.5 km 



Figure 4.8 Frequency Spectra at a depth of 150m and range of 2.5 Km 
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Figure 4. Depths 150 m. Ranges3.0 km 



Figure 4.9 Frequency Spectra at a depth of 150m and range of 3.0 Km 


X -j o'* Figure 5. Dep1hs200 m. Ran9e>2.5 km 



Figure 4.10 Frequency Spectra at a depth of 200m and range of 2.5 Km 
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,(in'* Figure 6. Dep<h=200 m. RangesJ.O km 



Figure 4.11 Frequency Spectra at a depth of 200m and range of 3.0 Km 


xio"* Figure 7. Depth*250 m. Range=2.5 km 



Figure 4.12 Frequency Spectra at a depth of 250m and range of 2.5 Km 
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Figure 8. Oep1h=250 m. Range-3.0 km 



Figure 4.13 Frequency Spectra at a depth of 250m and range of 3.0 Km 


xIO"® Figure 9. Deplh=300 m. Ran9e=2,5 km 



Figure 4.14 Frequency Spectra at a depth of 300m and range of 2.5 Km 
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K)-* Figure 10. Depth=300 m. Range=3.0 km 



Figure 4.15 Frequency Spectra at a depth of 300m and range of 3.0 Km 


Figure 11. Depth=3S0 m. Range=2.S km 



Figure 4.16 Frequency Spectra at a depth of 350m and range of 2.5 Km 
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xio"* Figure 12. Depth=350 m. Ran9e=3.0 km 



Figure 4.17 Frequency Spectra at a depth of 350m and range of 3.0 Km 

The main program, PEOUTMAIN was modified to compute average spectra for 
ranges of 2.5 and 3.0 km over depths of 100, 150, 200, 150, 300, and 350 m. Figures 4.18 
through 4.25 are the results. Figures 4.24 and 4.25 show spectra averaged over the six 
depths at ranges of 3.5 and 4.0 km respectively. 
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Figure 4.18 Average of Spectra at Ranges 2.5 and 3.0 km for Depth of 100 m 
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Figure 4.21 Average of Spectra at Ranges 2.5 and 3.0 km for Depth of 250 m 
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Figure 4.22 Average of Spectra at Ranges 2.5 and 3.0 km for Depth of 300 m 
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Figure 4.24 Average of Spectra at Depths 100-350 m at Range of 2.5 km 



Figure 4.25 Average of Spectra at Depths 100-350 m at Range of 3.0 km 
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V. ADVANCED MODELING ANALYSIS 


In the previous chapter we examined the influence from a dynamically moving 
surface and verified that the updated model performs adequately. We now wish to 
examine the performance of the code in the presence of more realistic, complicated 
surfaces. To do that we shall consider a rough surface based on an empirical fetch-limited 
ocean wave spectrum and how such modeling compares with measured data. 

A. PIERSON-MOSKOWITZ SPECTRUM 

Pierson and Moskowitz (1964) analyzed wave observations in the North Atlantic 
and proposed a form of energy spectrum of fully developed waves for each wind speed. 
They assumed that if the wind blew steadily for a sufficiently long time and over a 
sufficiently large area, the waves would come into equilibrium. They calculated the 
surface wave spectra for various wind speeds and determined that the spectra were of the 
form 


2 



4 

S(o)) = —^exp 

-P 



a> 


K^(0 ) 



(5.1) 


where co = , f is the surface wave frequency in Hertz, a = %.lx 10^^, P = 0.74, and 

cOq = g IU , where U is the wind speed reported by weather ships at a height of 19.5 

meters above the sea surface. Equation (5.1) is shown plotted in Fig. 5.1 with wind 
speeds from 15 to 20 m/s. Note the shift in peak response to slightly higher wave 
frequencies as the winds diminish. 
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B. JONSWAP SPECTRUM 

Hasselmann et al. (1973) determined after analyzing data collected during the 
Joint North Sea Wave Observation Project (JONSWAP) conducted in 1973 that the wave 
spectrum is never fully developed. It was found that for all wind speeds, the energy 
spectrum for growing waves took on similar characteristics depending on the stage of 
development of the waves and that this was due to the effect of non-linear wave 
interactions. A wave spectrum was proposed of the following form 


S(co) = ag^co ^ exp 

where S is the peak enhancement factor, 



r ^ 

-4 


CO 


-1.25 




CO 



\ p J 





(5.2) 


S = exp 



2(7q CO 


2 

p 


(5.3) 
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with y = 3.3, = 0.07 for co < co^, and = 0.09 for co> co^, while a is a function 

of fetch, X, and wind speed, U. Fetch is the distance from the lee shore or the distance 
over which the wind blows with a constant velocity, such that 


a = 0.076 


V U j 


The peak frequency is calculated by 




0 )^ = ln 


j\U‘ j 


(5.4) 


(5.5) 


Equation (5.2) is plotted in Fig. 5.2 with a constant fetch of 7.5 km and wind 
speeds varying from 10 to 30 m/s. Note how the peak frequencies diminish in amplitude 
as the wind speed diminishes and how the diminishing peaks increase in frequency. 
Figure 5.3 is a plot of (5.2) with wind speed constant at 10 m/s and the fetch varying from 
4 to 12 km. Note further how the amplitude of the peak frequencies increases with 
increasing fetch. Thus one can conclude that in coastal regions, the wind acts on a limited 
fetch with the result that the sea swell components may be less than what otherwise might 
occur in an open ocean environment and occur at higher frequency. 
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Figure 5.2 JONSWAP Spectrum with Varying Wind Speed and Constant Fetch 



Figure 5.3 JONSWAP Spectrum with Varying Fetch and Constant Wind Speed 
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c. 


HFA97 EXPERIMENT 


In September 1997, a high frequency acoustic experiment (HFA97) was 
conducted in a central region of the Delaware Bay. The intent of the HFA97 experiment 
was to measure cause and effect between the ocean environment and acoustic 
propagation. Two stable tripods, each having an acoustic source and three receiving 
hydrophones, were placed in 15 m of water separated by 387 m. On each tripod the 
source was located 3.125 m above the sea floor and the three receiving hydrophones were 
located at 0.33, 1.33 and 2.18 m. See Fig. 5.4 for details about the experimental setup. 

A broadband chirp signal in the frequency range from 1-18 kHz was transmitted 
every 0.345 seconds for an interval in one case of 5 seconds repeated every 10 minutes 
and in a second case for an interval of 40 seconds repeated every hour. 

In prior analysis of HFA97 data, remotely received signals across the three 
hydrophones were used with a beamforming technique to calculate signal arrival angle as 
a function of arrival time (Badiey, et al, 2000). By considering the geometry of the 
HFA97 experimental setup, the resulting beamformed plots can be used to easily 
distinguish that portion of the received signal corresponding to single surface reflected 
wave paths, as shown in Fig. 5.5. 
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Experimental Setup 
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Figure 5.4 HFA97 Experimental Setup, numbers 2, 3, 4, and 5 represent single 
surface reflected ray paths (After Ref. [Heitsenrether, Badiey, 2004]). 
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Figure 5.5 Received signal arrival angle versus arrival time for a calm 
period. Single surface reflected ray paths are easily distinguished in the 
signal (From Ref. [Heitsenrether, Badiey, 2004]). 
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D. IMPLEMENTATION OF JONSWAP INTO MMPE MODEL 


A sea surface deviation wavenumber spectrum, W(k) obtained from the 
JONSWAP wave spectrum, Equation (5.2), has the form 


W(k) = 


do) 

dk 


(5.6) 


With the wave dispersion relationship a) = 


^|kg , 


dcD _ 1 ^ 

dk 2 ]l k 


(5.7) 


where k is the wavenumber of ocean waves. From Equations (5.2) through (5.7), the sea 
surface wavenumber spectrum can be shown to have the form 


W(k) = 



f otg^ikg) 

k 


1 


exp 


ak 


-3 


(Op\kgy 


(3.3)" 


(5.8) 


From the wavenumber spectrum (5.8), code was developed for implementation 
into the rough surface forward scatter subroutine, ZSGEN, of the MMPE. Recall from 
Section III that ZSGEN calculates a roughness spectrum that serves as an input to 
ENVPROP for calculating the image ocean potential function. The changes to the 
ZSGEN program include the swapping of the one dimensional rough surface wave 
number spectrum calculation with the JONSWAP sea surface wavenumber spectrum 
calculation (5.8). The input parameters were selected from the setup in Figure 5.4 and 
include a source depth of 12 m, a maximum range of 0.387 km, and a water depth of 15 
m. An isovelocity sound speed of 1500 m/s was selected, with sea bottom sound speed of 

1650 m/s, and sea bottom density of 1.9 g/cc. A relatively high sea bottom attenuation of 
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0.5 db/km/Hz was chosen so as reduce scatter from the sea floor. A fetch distance of 7.5 
km was input to the model as an estimate of actual conditions during HFA97. The sound 
source inputs included a center frequency of 9 kHz, a bandwidth of 16 kHz, and 256 
frequencies selected for processing. 

The model was run for various wind speeds. Figures 5.6 through 5.9 are the 
results for wind speeds of 1, 5, 10, and 15 m/s. When comparing these figures with Fig. 
5.5 the similarities between the model and the empirical data are apparent. Note the 
dispersion of energy with the surface reflected paths as the wind speed increases and how 
the direct path is unaffected. 



Time (sec) 


Figure 5.6 MMPE results of arrival angle versus arrival time for a 

wind speed of 1 m/s. 
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Figure 5.7 MMPE results of arrival angle versus arrival time for a 

wind speed of 5 m/s 
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Figure 5.8 MMPE results of arrival angle versus arrival time for a 

wind speed of 10 m/s. 
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Figure 5.9 MMPE results of arrival angle versus arrival time for a 

wind speed of 15 m/s. 


E. COMPUTING STANDARD DEVIATION OF ARRIVAL TIME AND 
ARRIVAL ANGLE 


Two MATLAB programs were written for calculating the standard deviation of 
arrival angle and arrival time (see Appendix A for details of MATLAB code). For both 
algorithms, a narrow window or slice is sampled from the single surface return 
transmission loss pressure data for positive angles from 0 to 39.5 degrees and time width 
of 0.36 msecs. The minimum transmission loss value within the window is found and 
established as a reference point. Next, the data is sorted into twelve bins of transmission 
loss data each of 1 dB difference referenced to the minimum transmission loss value. 
The data in the 12 bins are next weighted by values of 32, 25, 20, 16, 13, 10, 8, 6, 5, 4, 3, 
and 3. An emphasis weighting value is added to the angle of arrival algorithm and 
operates on the higher transmission loss values. The emphasis adds linearity at the higher 
wind speeds 
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Measuring dispersion within a window or slice that contains data represented by 
position and magnitude, both of which are changing (dispersing) as the wind speed 
increases, is a challenging task. Various combinations involving window sizes, number 
of thresholds, threshold levels, and various weighting associated with each threshold 
were examined while trying to maintain a rational basis for the algorithms. Although not 
exact, the algorithms are considered a good approximation with the results shown in Figs. 
5.10 and 5.11 being optimum when compared to HFA97 measured data results shown in 
Figs. 5.12 and 5.13. 
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• MMPE Model Data 
- Linear Regression 


Figure 5.10 Standard Deviation of Arrival Angle vs. Wind Speed for MMPE Model 

Data 
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Figure 5.11 Standard Deviation of Arrival Time vs. Wind Speed for MMPE Model 

Data 



Figure 5.12 Standard Deviation of Arrival Angle vs. Wind Speed for Measured 
HFA97 data (After Ref. [Heitsenrether, Badiey, 2004]). 
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Figure 5.13 Standard Deviation of Arrival Time vs. Wind Speed for Measured 
HFA97 data (After Ref. [Heitsenrether, Badiey, 2004]). 

Figures 5.14 and 5.15 compare regression analysis of measured data with 
modeled data for both arrival angle and arrival time (see Appendix B for details about 
regression statistics results). What can be gathered from the regression analysis is not 
only that acoustical energy dispersion increases with increasing wind speed, but also how 
the steepness of the slope of measured data and modeled data relates to the rate of 
acoustical energy dispersion. Thus, the MMPE model is predicting slightly lower energy 
dispersion rates when compared to an analysis of measured data. Note, however, that the 
MATLAB algorithms that compute the MMPE data arrival angle and arrival time 
statistics are considered approximations. Also, the window or slice of data that the 
algorithm examines may include a small portion of other acoustic energy besides 
reflection from just a single surface bounce. 

Another model employed as a combination of an empirical sea surface model and 
a Gaussian beam tracing model (BELLHOP), termed the BELLHOP/JONSWAP model, 
was compared with measured data from the HPA97 experiment (Heitsenrether, Badiey, 
2004). Eigures 5.16 and 5.17 show results of the BELLHOP/JONSWAP model with 
input parameters taken from the HPA97 experiment. 
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Figure 5.14 Regression Analysis Comparing Standard Deviation of Arrival Angle vs. 
Wind Speed for Measured and Modeled Data 



Figure 5.15 Regression Analysis Comparing Standard Deviation of Arrival Time vs. 
Wind Speed for Measured and Modeled Data 
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Figure 5.16 Standard Deviation of Arrival Angle vs. Wind Speed for 
BELLHOP/JONSWAP Model data (After Ref. [Heitsenrether, Badiey, 2004]). 



Figure 5.17 Standard Deviation of Arrival Time vs. Wind Speed for 
BELLHOP/JONSWAP Model data (After Ref. [Heitsenrether, Badiey, 2004]). 
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Finally, Figures 5.18 and 5.19 show a regression analysis comparison of the 
BELLHOP/JONSWAP model, measured data from the HFA97 experiment, and the 
MMPE model for both the standard deviation of arrival angle vs. wind speed, and the 
standard deviation of arrival time vs. wind speed. By examining the slopes of the three 
curves it is interesting to note how the BELLHOP/JONSWAP model is predicting a 
slightly greater dispersion rate when compared to measured data and how the MMPE 
model is predicting a slightly less dispersion rate than measured data. Eurthermore, the 
BELLHOP/JONSWAP model has a negative intercept with the y-axis, whereas the 
MMPE has a positive intercept, and the measured data extrapolates to approximately 
zero. 


Regression Comparison with BELLOP/JONSWAP, Measured, and 
MMPE Modeled Data for Stdv of Arrival Angle vs. Wind Speed 
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Figure 5.18 Regression Analysis Comparing Standard Deviation of Arrival Angle vs. 
Wind Speed for BELLHOP/JONSWAP, Measured, and MMPE Modeled Data 
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Regression Comparison with BELLHOP/JONSWAP, Measured, and 
MMPE Modeled Data for Stdv Arrival Time vs Wind Speed 
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Figure 5.19 Regression Analysis Comparing Standard Deviation of Arrival Time vs. 

Wind Speed for BELLHOP/JONSWAP, Measured, and MMPE Modeled Data 

Although the details of the BELLHOP/JONSWAP model and corresponding data 
analysis are not available here, we can make several comments on the issues related to 
the MMPE implementation. It was previously noted that the MMPE surface scatter 
implementation requires an ad-hoc minimum surface displacement sampling of 1/5. The 
effect this has on the results is unknown, but may be expected to introduce the most 
errors at small surface displacements (i.e., small wind speeds). It was also noted that the 
variability analysis may capture numerical noise which could also affect the low wind 
speed results. Both of these could contribute to the non-zero variability intercept of the 
MMPE results. 

In spite of these potential issues, the MMPE/JONS WAP model produced results 
consistent with the measured data and, arguably, as good as the BELLHOP/JONSWAP 
model predictions. It is also useful to note some comparison between the BELLHOP and 
MMPE predictions. In both cases, the predicted rate of increase of variability in both 
arrival angle and arrival time is higher using the BELLHOP model than the MMPE 
model. This could be due to the higher sensitivity of the ray-based model, although there 
is no direct evidence of this. Eurthermore, the predicted rate of increase of variability in 



BELLHOP/JONSWAP 

Model 

Measured Data 
MMPE Modeled Data 


5.00 10.00 15.00 

Wind Speed (m/s) 


20.00 
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the arrival time from the MMPE model matched the measured data quite well, while the 
slope of the BELLHOP results is significantly larger. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


The littoral waters are challenging acoustically with sea surface motion constantly 
changing. Developing models that can predict the physical acoustic behavior of this 
medium with some significant degree of accuracy is a complex effort. The approach 
taken by this thesis is to build upon existing work related to the Monterey Miami 
Parabolic Equation model. In each of the chapters in this thesis, analysis has been 
presented followed by model implementation discussions and results. 

Background information about the MMPE was provided as a means to establish a 
baseline upon which the advanced work of developing sea surface models could be built. 
A static rough surface forward scatter model was developed which displayed increases in 
sea surface acoustic scatter as the rms roughness input parameter was increased. A CW 
test of scattering from a sinusoidal surface was conducted with the results showing 
reflected rays due to Bragg scatter traveling in distinct directions, as expected. As the 
surface was allowed to move dynamically, the directions or angles changed but the 
distinction remained evident. Analysis was conducted on the scattered field from the 
moving sinusoidal surface and Doppler shifts were determined to be in agreement with 
calculated Bragg scatter conditions. 

In Section V an empirical fetch-limited ocean wave spectrum (JONSWAP) was 
presented. Model parameters from the HPA97 experiment were used as input and 
compared with measured data. The results look promising and clearly show dispersion of 
energy associated with single surface reflected rays as wind speed increases, and clearly 
show how the direct path (non-surface interaction) remains unchanged. A comparison of 
linear regression analysis of measured data with MMPE modeled data shows that not 
only acoustical energy dispersion increases with increasing wind speed, but also how the 
steepness of the slope of the measured data and modeled data relates to the rate of 
acoustical energy dispersion. The MMPE/JONSWAP predictions of energy dispersion 
rates seemed to compare quite favorably for the variability of acoustic travel time. 
However, the predicted rate for arrival angle variability was too low. 
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The two MATLAB algorithms that compute arrival angle and arrival time 
statistics for the MMPE model data are considered good approximations but are not 
exact. Furthermore, the MMPE results had higher accuracy in arrival time measurements 
than arrival angle. A combination of these issues, and perhaps others, may be the cause of 
the underestimate of energy dispersion in predicted arrival angle. 

Finally, a regression analysis comparison was made of the BELLHOP/JONSWAP 
model, the measured data from the HFA97 experiment, and the MMPE model for both 
the standard deviation of arrival angle and the standard deviation of arrival time vs. wind 
speed. By examining the slopes of the three curves, the BELLHOP/JONSWAP model 
predicts a greater dispersion rate than the MMPE model in both cases. Furthermore, the 
BELLHOP/JONSWAP model was observed to have a negative intercept with the y axis, 
whereas the MMPE had a positive intercept and the measured data extrapolates to 
approximately zero. It is unclear why both numerical models missed a physically 
plausible zero intercept. 

The results presented show that surface scatter can significantly affect arrival 
time, arrival angle, and frequency spread that leads to degradation of acoustic signal 
coherence. Acoustical communication systems can benefit from models predicting the 
environmental conditions affecting coherence and potentially these systems can adjust or 
compensate for the existing conditions. 

Some recommendations for future work include examining Doppler shifts for the 
dynamically evolving rough surface, improving upon the window/slice statistical 
MATLAB algorithms or examining another approach for determining standard deviation 
of arrival angle and arrival time, and validating the algorithms with additional measured 
data. 
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APPENDIX A. DATA PROCESSING 


MATLAB Program for Computing Standard Deviation of Arrival Angle 

%CompAngleSTDVrev6st.m computes the standard deviation of arrival angle in 
degrees. 

%Runs in conjunction with PEOUT2 

clear tlpressbeamtstl; 
clear tlpressbeamtst2; 

for n6=138:1:140; % Corresponds to window slice of 0.00012 seconds per 

increment or 0.00036 secs. 

for n5=180:l:240; % by positive angle of zero to 39.5 degrees 
% Grab slice 
nl0=n5-179; 
nll=n6-137; 

tlpressbeamtst2(n 10,n 11 )=tlpressbeam(n5 ,n6); 
end 
end 

Amin=min(tlpressbeamtst2); % Find least value of TL 
Bmin=min(Amin) 

threshold 1 l=Bmin+l; % thresholds established for 12 bins of 1 db ref to min TL 

threshold 12=Bmin+2; 

threshold 13=Bmin+3; 

threshold 14=Bmin+4; 

threshold 15=Bmin+5; 

threshold 16=Bmin+6; 

threshold 17=Bmin+7; 

threshold 18=Bmin+8; 

threshold 19=Bmin+9; 

threshold20=Bmin+10; 

threshold21 =Bmin+11; 

threshold22=Bmin+12; 


emp=20; % emphasis value placed on higher TL values; Adds linearity to 

computation 

n7=0; % at higher windspeeds 
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for n6=138:l:140; % Corresponds to window slice of 0.00036 seconds 
for n5= 180:1:240; % by positive angle of zero to 39.5 degrees 

nl0=n5-179; 

nll=n6-137; 

tlpressbeamtst 1 (n 10,n 11 )=tlpressbeam(n5 ,n6); 
if tlpressbeam(n5,n6)<thresholdl 1; 

n8=n7+32; % weighting of 32 times for transmission loss level below Bmin 

+ 1 db 

if n7==0; 

n7=l; 

end 

for n9=n7:l:n8; 

Etest(n9)=theta(n5); 

end 

n7=n8; 

else 

if tlpres sbeam(n5 ,n6)<threshold 12; 
n8=n7+25; % weighting of 25 times 
if n7==0; 

n7=l; 

end 

for n9=n7:l:n8; 

Etest(n9)=theta(n5); 

end 

n7=n8; 

else 

if tlpres sbeam(n5 ,n6 )<threshold 13; 
n8=n7+20; % weighting of 20 
if n7==0; 

n7=l; 

end 

for n9=n7:l:n8; 

Etest(n9)=theta(n5); 

end 

n7=n8; 

else 

if tlpres sbeam(n5 ,n6 )<threshold 14; 
n8=n7+16; % weighting of 16 
if n7==0; 

n7=l; 

end 

for n9=n7:l:n8; 
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Etest(n9)=theta(n5); 

end 

n7=n8; 

else 

if tlpressbeam(n5,n6)<thresholdl5; 
n8=n7+13; % weighting of 13 
if n7==0; 

n7=l; 

end 

for n9=n7:l:n8; 

Etest(n9)=theta(n5); 

end 

n7=n8; 

else 

if tlpressbeam(n5,n6)<thresholdl6; 
n8=n7+10; % weighting of 10 
if n7==0; 

n7=l; 

end 

for n9=n7:l:n8; 

Etest(n9)=theta(n5); 

end 

n7=n8; 

else 

if tlpressbeam(n5,n6)<thresholdl7; 

n8=n7+8+emp; % weighting of 8 plus emphasis value 
if n7==0; 

n7=l; 

end 

for n9=n7:l:n8; 

Etest(n9)=theta(n5); 

end 

n7=n8; 

else 

if tlpressbeam(n5,n6)<thresholdl8; 

n8=n7+6+emp; % weighting of six plus emphasis value 
if n7==0; 

n7=l; 

end 

for n9=n7:l:n8; 

Ete st(n9)=theta(n5); 
end 
n7=n8; 
else 
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if tlpressbeam(n5 ,n6)<threshold 19; 

n8=n7+5+emp; % weighting of 5 plus emphasis value 
if n7==0; 

n7=l; 

end 

for n9=n7:l:n8; 

Etest(n9)=theta(n5); 

end 

n7=n8; 

else 

if tlpressbeam(n5,n6)<threshold20; 

n8=n7+4+emp; % weighting of 4 plus emphasis value 
if n7==0; 

n7=l; 

end 

for n9=n7:l:n8; 

Etest(n9)=theta(n5); 

end 

n7=n8; 

else 

if tlpressbeam(n5,n6)<threshold21; 

n8=n7+3+emp; % weighting of 3 plus emphasis 

value 

if n7==0; 

n7=l; 

end 

for n9=n7:l:n8; 

Etest(n9)=theta(n5); 

end 

n7=n8; 

else 

if tlpressbeam(n5,n6)<threshold22; 

n8=n7+3+emp; % weighting of 3 plus emphasis 

value 

if n7==0; 

n7=l; 

end 

for n9=n7:l:n8; 

Etest(n9)=theta(n5); 

end 

n7=n8; 

end 

end 

end 
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end 


end 

end 

end 

end 

end 


end 

end 

end 

end 

end 


Angtest=std(Etest)% compute standard deviation of arrival angle 
clear Etest 

MATLAB Program for Computing Standard Deviation of Time of Arrival 

%CompTimeSTDV computes the standard deviation of time of arrival and runs 
% in conjunction with PEOUT2 

clear tlpressbeamtstl; 
clear tlpressbeamtst2; 

for n6=138:1:140; % Corresponds to window slice of 0.00012 seconds per 

increment 

% or 0.00036 secs. 

for n5=180:l:240; % by positive angle of zero to 39.5 degrees 

nl0=n5-179; 

nll=n6-137; 

tlpressbeamtst2(n 10,n 11 )=tlpressbeam(n5 ,n6); 
end 
end 

Amin=min(tlpressbeamtst2); % Eind least value of TL 
Bmin=min(Amin) 

threshold 1 l=Bmin+l; % thresholds established for 12 bins of 1 db ref to min TL 
threshold 12=Bmin+2; 
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threshold 13=Bmin+3; 
threshold 14=Bmin+4; 
threshold 15=Bmin+5; 
threshold 16=Bmin+6; 
threshold 17=Bmin+7; 
threshold 18=Bmin+8; 
threshold 19=Bmin+9; 
threshold20=Bmin+10; 
threshold21 =Bmin+11; 
threshold22=Bmin+12; 


emp=0 % emphasis value, zero means no emphasis added 
n7=0; 


for n6=138:l:140; % Corresponds to window slice of 0.00036 seconds 
for n5= 180:1:240; % by positive angle of zero to 39.5 degrees 

nl0=n5-179; 

nll=n6-137; 

tlpressbeamtst 1 (n 10,n 11 )=tlpressbeam(n5 ,n6); 
if tlpressbeam(n5,n6)<thresholdl 1; 

n8=n7+32; % weighting of 32 times for transmission loss level below Bmin 

+ 1 db 

if n7==0; 

n7=l; 

end 

for n9=n7:l:n8; 

Ftest(n9)=timeout(n5); 

end 

n7=n8; 

else 

if tlpres sbeam(n5 ,n6)<threshold 12; 
n8=n7+25; % weighting of 25 times 
if n7==0; 

n7=l; 

end 

for n9=n7:l:n8; 

Ftest(n9)=timeout(n5); 

end 

n7=n8; 

else 
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if tlpres sbeam(n5 ,n6 )<threshold 13; 
n8=n7+20; % weighting of 20 times 
if n7==0; 

n7=l; 

end 

for n9=n7:l:n8; 

Ftest(n9)=timeout(n5); 

end 

n7=n8; 

else 

if tlpressbeam(n5 ,n6)<threshold 14; 
n8=n7+16; % weighting of 16 times 
if n7==0; 

n7=l; 

end 

for n9=n7:l:n8; 

Ftest(n9)=timeout(n5); 

end 

n7=n8; 

else 

if tlpressbeam(n5,n6)<thresholdl5; 
n8=n7+13; % weighting of 13 times 
if n7==0; 

n7=l; 

end 

for n9=n7:l:n8; 

Ftest(n9)=timeout(n5); 

end 

n7=n8; 

else 

if tlpressbeam(n5,n6)<thresholdl6; 
n8=n7+10; % weighting of 10 times 
if n7==0; 

n7=l; 

end 

for n9=n7:l:n8; 

Ftest(n9)=timeout(n5); 

end 

n7=n8; 

else 

if tlpressbeam(n5,n6)<thresholdl7; 
n8=n7+8+emp; % weighting of 8 times 
if n7==0; 

n7=l; 

end 
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for n9=n7:l:n8; 

Ftest(n9)=timeout(n5); 

end 

n7=n8; 

else 

if tlpressbeam(n5,n6)<thresholdl8; 

n8=n7+6+emp; % weighting of 6 times 
if n7==0; 

n7=l; 

end 

for n9=n7:l:n8; 

Ftest(n9)=timeout(n5); 

end 

n7=n8; 

else 

if tlpressbeam(n5 ,n6)<threshold 19; 

n8=n7+5+emp; % weighting of 5 times 
if n7==0; 

n7=l; 

end 

for n9=n7:l:n8; 

Ftest(n9)=timeout(n5); 

end 

n7=n8; 

else 

if tlpressbeam(n5,n6)<threshold20; 

n8=n7+4+emp; % weighting of 4 times 
if n7==0; 

n7=l; 

end 

for n9=n7:l:n8; 

Ftest(n9)=timeout(n5); 

end 

n7=n8; 

else 

if tlpressbeam(n5,n6)<threshold21; 

n8=n7+3+emp; % weighting of 3 times 
if n7==0; 

n7=l; 

end 

for n9=n7:l:n8; 

Ftest(n9)=timeout(n5); 

end 

n7=n8; 

else 
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if tlpressbeam(n5,n6)<threshold22; 

n8=n7+3+emp; % weighting of 3 times 
if n7==0; 

n7=l; 

end 

for n9=n7:l:n8; 

Ftest(n9)=timeout(n5); 

end 

n7=n8; 

end 

end 

end 

end 

end 

end 

end 

end 

end 

end 

end 

end 

end 

end 


Timtest=std(Ftest)% compute standard deviation of arrival time 
clear Ftest 
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APPENDIX B. LINEAR REGRESSION STATISICAL ANALYSIS 


Standard Deviation of Arrival Angle vs. Wind Speed for MMPE Model Data 

SUMMARY OUTPUT 


Regression Statistics 


Multiple R 

0.986387 

R Square 

0.972959 

Adjusted R 


Square 

0.968452 

Standard Error 

0.093242 

Observations 

8 


ANOVA 



df 

SS 

MS 

F 

Significance 

F 

Regression 

1 

1.876923 

1.876923 

215.8833 

6.24E-06 

Residual 

6 

0.052165 

0.008694 



Total 

7 

1.929088 






Coefficients 

Standard 

Error 

tStat 

P-vaiue 

Lower 95% 

Upper 

95% 

Intercept 

0.450394 

0.066487 

6.77415 

0.000506 

0.287706 

0.613083 

X Variable 1 

0.102831 

0.006999 

14.69297 

6.24E-06 

0.085706 

0.119956 


RESIDUAL OUTPUT 


Observation 

Predicted Y 

Residuais 

1 

0.553225 

0.006775 

2 

0.758887 

-0.00889 

3 

0.964549 

0.035451 

4 

1.170211 

-0.13021 

5 

1.478704 

0.091296 

6 

1.684366 

0.005634 

7 

1.787197 

0.112803 

8 

1.992859 

-0.11286 
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Standard Deviation of Arrival Time vs. Wind Speed for MMPE Model Data 


SUMMARY OUTPUT 


Regression Statistics 


Multiple R 

0.977426 

R Square 

0.955361 

Adjusted R 


Square 

0.947921 

Standard Error 

0.106244 

Observations 

8 


ANOVA 



df 

SS 

MS 

F 

Significance 

F 

Regression 

1 

1.449474 

1.449474 

128.4118 

2.83E-05 

Residual 

6 

0.067726 

0.011288 



Total 

7 

1.5172 






Coefficients 

Standard 

Error 

tStat 

P-vaiue 

Lower 95% 

Upper 

95% 

Intercept 

0.464479 

0.075758 

6.131106 

0.000861 

0.279106 

0.649852 

X Variable 1 

0.090366 

0.007975 

11.33189 

2.83E-05 

0.070853 

0.109879 


RESIDUAL OUTPUT 


Observation 

Predicted Y 

Residuais 

1 

0.554845 

-0.02485 

2 

0.735577 

-0.09558 

3 

0.91631 

0.10369 

4 

1.097042 

-0.05704 

5 

1.368141 

0.161859 

6 

1.548873 

0.021127 

7 

1.639239 

0.020761 

8 

1.819972 

-0.12997 
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Standard Deviation of Arrival Angle vs. Wind Speed for Measured HFA97 Data 

SUMMARY OUTPUT 


Regression Statistics 


Multiple R 

0.832309 

R Square 

0.692738 

Adjusted R 


Square 

0.679378 

Standard Error 

0.439401 

Observations 

25 


ANOVA 



df 

SS 

MS 

F 

Significance 

F 

Regression 

1 

10.01175 

10.01175 

51.85459 

2.48E-07 

Residual 

23 

4.44069 

0.193073 



Total 

24 

14.45244 






Coefficients 

Standard 

Error 

tStat 

P-vaiue 

Lower 95% 

Upper 

95% 

Intercept 

0.155232 

0.21186 

0.732711 

0.471137 

-0.28303 

0.593498 

X Variable 1 

0.171888 

0.02387 

7.201013 

2.48E-07 

0.122509 

0.221267 


RESIDUAL OUTPUT 


Observation 

Predicted Y 

Residuais 

1 

0.385562 

-0.38556 

2 

0.385562 

-0.18556 

3 

0.316807 

-0.11681 

4 

1.038737 

-0.09874 

5 

0.923572 

-0.09357 

6 

0.923572 

0.206428 

7 

0.923572 

0.386428 

8 

1.153902 

0.346098 

9 

1.315476 

0.184524 

10 

1.461581 

0.413419 

11 

1.499397 

-0.5594 

12 

1.499397 

-0.3694 

13 

1.691911 

-0.19191 

14 

1.769261 

0.480739 

15 

1.844892 

0.235108 

16 

1.844892 

-0.71489 

17 

1.922241 

0.327759 

18 

1.922241 

0.707759 

19 

2.07694 

-0.19694 
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20 

2.152571 

0.097429 

21 

2.152571 

0.097429 

22 

2.152571 

0.657429 

23 

2.30727 

-0.05727 

24 

2.460251 

-1.15025 

25 

2.460251 

-0.02025 


Standard Deviation of Arrival Time vs. Wind Speed for Measured HFA97 Data 

SUMMARY OUTPUT 


Regression Statistics 


Multiple R 

0.877373 

R Square 

0.769784 

Adjusted R 


Square 

0.757668 

Standard Error 

0.217574 

Observations 

21 


ANOVA 



df 

SS 

MS 

F 

Significance 

F 

Regression 

1 

3.007481 

3.007481 

63.53124 

1.77E-07 

Residual 

19 

0.899434 

0.047339 



Total 

20 

3.906914 






Coefficients 

Standard 

Error 

tStat 

P-vaiue 

Lower 95% 

Upper 

95% 

Intercept 

-0.02323 

0.114162 

-0.20351 

0.840896 

-0.26218 

0.21571 

X Variable 1 

0.108426 

0.013603 

7.970649 

1.77E-07 

0.079954 

0.136898 


RESIDUAL OUTPUT 


Observation Predicted Y Residuais 


1 

0.098204 

0.181796 

2 

0.145911 

0.134089 

3 

0.170849 

-0.17085 

4 

0.461432 

-0.27143 

5 

0.461432 

-0.18143 

6 

0.485285 

0.174715 

7 

0.606723 

-0.13672 

8 

0.704306 

-0.14431 
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9 

0.800806 

-0.05081 

10 

0.800806 

-0.05081 

11 

0.849598 

0.090402 

12 

0.946097 

-0.0061 

13 

0.994889 

0.225111 

14 

1.018742 

-0.17874 

15 

1.067534 

-0.03753 

16 

1.067534 

0.242466 

17 

1.04303 

0.64697 

18 

1.188972 

-0.05897 

19 

1.212825 

-0.08283 

20 

1.334263 

-0.11426 

21 

1.430762 

-0.21076 




Standard Deviation of Arrival Angle vs. Wind Speed for BELLHOP/JONSWAP 
Model Data 


Multiple R 

0.983743 

R Square 

0.96775 

Adjusted R 


Square 

0.965269 

Standard Error 

0.169868 

Observations 

15 


ANOVA 



df 

SS 

MS 

F 

Significance 

F 

Regression 

1 

11.25626 

11.25626 

390.0964 

4.45E-11 

Residual 

13 

0.375116 

0.028855 



Total 

14 

11.63137 






Coefficients 

Standard 

Error 

tStat 

P-vaiue 

Lower 95% 

Upper 

95% 

Intercept 

-0.33638 

0.127809 

-2.63193 

0.020711 

-0.6125 

-0.06027 

X Variable 1 

0.271122 

0.013727 

19.75086 

4.45E-11 

0.241466 

0.300778 


RESIDUAL OUTPUT 


Observation Predicted Y Residuals 
1 0.476981 -0.19698 
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2 

0.699301 

-0.1393 

3 

1.14123 

-0.02123 

4 

1.439464 

-0.03946 

5 

1.58587 

0.19413 

6 

1.808189 

0.001811 

7 

1.881392 

0.068608 

8 

2.030509 

0.059491 

9 

2.103712 

0.126288 

10 

2.326032 

0.043968 

11 

2.472438 

0.037562 

12 

2.767961 

0.162039 

13 

2.917078 

0.152922 

14 

3.212601 

-0.0026 

15 

3.657241 

-0.44724 


Standard Deviation of Arrival Time vs. Wind Speed for BELLHOP/JONSWAP 
Model Data 

Standard 

Error 0.16319 

Observations 18 


ANOVA 



df 

SS 

MS 

F 

Significance 

F 

Regression 

1 

7.706266 

7.706266 

289.3727 

1.14E-11 

Residuai 

16 

0.426095 

0.026631 



Totai 

17 

8.132361 






Coefficients 

Standard 

Error 

tStat 

P-vaiue 

Lower 95% 

Upper 

95% 

Intercept 

-0.39017 

0.100389 

-3.88659 

0.00131 

-0.60299 

0.17736 

X Variable 1 

0.179012 

0.010523 

17.01096 

1.14E-11 

0.156704 

0.20132 


RESIDUAL OUTPUT 


Observation Predicted Y Residuals 

1 0.205938 -0.11594 

2 0.256062 -0.06606 

3 0.306185 -0.11618 

4 0.404641 -0.03464 
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5 

0.605135 

-0.04513 

6 

0.703592 

0.036408 

7 

0.902295 

0.117705 

8 

1.002542 

0.117458 

9 

1.201245 

0.008755 

10 

1.299701 

0.000299 

11 

1.399948 

0.090052 

12 

1.500195 

-0.01019 

13 

1.598651 

0.071349 

14 

1.698898 

0.341102 

15 

1.897601 

0.142399 

16 

1.996058 

0.043942 

17 

2.096305 

-0.1463 

18 

2.295008 

-0.43501 
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